#
# Copyright (c) 2025 CAX Conversion Project
#

import numpy as np
from scipy.fft import fft
import numpy.fft as nf
from scipy.stats.distributions import chi2
import pywt
from datetime import datetime

def myfft(info):
    # information = json.loads(info)
    # print(information)
    time = info["time"]
    x = info["value"]
    fft_x = fft(x)                                            #  fft计算
    amp_x = abs(fft_x)/len(x)*2                                 # 纵坐标变换
    label_x = np.linspace(0, int(len(x)/2)-1, int(len(x)/2))    # 生成频率坐标
    amp = amp_x[0:int(len(x)/2)]                              # 选取前半段计算结果即可
    # amp[0] = 0                                              # 可选择是否去除直流量信号
    t1 = datetime.strptime(time[2], "%H:%M:%S:%f")
    t2 = datetime.strptime(time[1], "%H:%M:%S:%f")
    t = t1 - t2
    print(t)
    hh, mm, ss = map(str, str(t).split(':'))
    ss, fs = map(int, str(ss).split('.'))
    t = 1000*(ss + 60*(int(mm) + 60*int(hh)))+fs
    fs = 1 / t                                     # 计算采样频率
    fre = label_x/len(x)*fs                                   # 频率坐标变换
    pha = np.unwrap(np.angle(fft_x))                          # 计算相位角并去除2pi跃变
    # return amp, fre
    return amp.tolist()


def fft(data):
    sample_rate = int(data["frequency"])
    data = data["value"]
    data = list(map(float, data))
    times = np.arange(len(data)) / len(data)
    freqs = nf.fftfreq(times.size, times[1] - times[2])
    complex_array = nf.fft(data)
    pows = np.abs(complex_array)
    x = (freqs[freqs > 0])[:int(sample_rate / 2)]
    y = (pows[freqs > 0])[:int(sample_rate / 2)]
    x = np.array(x, dtype=int)
    y = np.round(y, 3)
    f_data = y
    dataDict = {"t": x.tolist(), "amp": f_data.tolist()}
    return dataDict


def bo(data):
    sampling_rate = 1024
    wavename = 'cgau8'
    totalscal = 256
    fc = pywt.central_frequency(wavename)
    cparam = 2 * fc * totalscal
    scales = cparam / np.arange(totalscal, 1, -1)
    [cwtmatr, frequencies] = pywt.cwt(data, scales, wavename, 1.0 / sampling_rate)
    dataDict = {"amp": cwtmatr.tolist(), "fre": frequencies.tolist()}
    return dataDict


def cwt(info):

    x = info["value"]
    x = list(map(float, x))
    wavelet = 'cgau8'
    w = 5
    z = 30
    fs = int(info["frequency"])
    time = 2
    f = w * z
    t = np.linspace(0, time - 1 / fs, int(time * fs))
    # x = (1 + 1 * np.sin(2 * np.pi * 20 * t)) * np.sin(2 * np.pi * f * t)
    # amp, fre = cwt(info, fs, 512, 'morl')
    totalscal = 1

    if wavelet not in pywt.wavelist():
        print('小波函数名错误')
    else:
        wfc = pywt.central_frequency(wavelet=wavelet)
        a = 2 * wfc * totalscal / (np.arange(totalscal, 0, -1))
        period = 1.0 / fs
        [cwtmar, fre] = pywt.cwt(x, a, wavelet, period)
        amp = abs(cwtmar)

    dataDict = {"t": t.tolist(), "amp": amp[0].tolist(), "fre": fre.tolist()}
    return dataDict

'''
功率谱分析
输入：
x：需要分析的时间序列(原始序列，未标准化或距平处理)
m：最大滞后相关长度，m取值范围最好在(n/10)~(n/3)之间，n为样本数，可以多次调整m获得最佳效果，通常取m=n/3
alpha1：红噪音检验信度
alpha2：白噪音检验信度
输出：
l：功率谱图的X坐标，对应的周期为2m/l，使用时自己调整tick labels
Sl：功率谱估计值
Sr：红噪音
Sw：白噪音
r1：落后一个时刻的自相关函数，用于查看使用哪种噪音检验
'''
def specx(data):
    data = data["value"]
    data = list(map(float, data))
    x = np.array(data)
    m = x.shape[0] // 3
    alpha1 = 0.05
    alpha2 = 0.05
    n = x.shape[0] # 序列的长度
    if np.std(x) != 0.0:
        x = (x - np.mean(x)) / np.std(x) # 标准化
    r1 = np.zeros((n-6))
    r2 = np.zeros((n-7))
    for i in np.arange(0, n-6):
        r1[i] = np.sum(x[:n-i]*x[i:])/x[:n-i].shape[0]
    for i in np.arange(1, n-6):
        r2[i-1] = np.sum(x[:n-i]*x[i:])/x[:n-i].shape[0]
    r2 = r2[::-1]
    r = np.hstack((r2, r1))
    l = np.arange(0, m+1, 1)
    tao = np.arange(1, m, 1)
    Sl  = np.zeros((m+1))
    Tl  = np.zeros((m+1))
    S0l = np.zeros((m+1))
    a = np.array((r.shape[0]+1)/2).astype('int32')
    r = r[a-1:a+m]
    a = r[1:-1]*(1+np.cos(np.pi*tao/m))
    for i in np.arange(2, m+1, 1):
        Sl[i-1] = (r[0]+np.sum(a*np.cos(l[i-1]*np.pi*tao/m)))/m
    Sl[0] = (r[0]+np.sum(a*np.cos(l[0]*np.pi*tao/m)))/(2*m)
    Sl[-1] = (r[0]+np.sum(a*np.cos(l[-1]*np.pi*tao/m)))/(2*m)
    for i in range(l.shape[0]):
        Tl[i] = 2*m/l[i]
    f = (2*n-m/2)/m
    S = np.mean(Sl)
    for i in range(l.shape[0]):
        S0l[i] = S*(1-r[1]*r[1])/(1+r[1]*r[1]-2*r[1]*np.cos(l[i]*np.pi/m))
    x2r = chi2.ppf(1-alpha1, df=f)
    Sr = S0l*x2r/f
    x2w = chi2.ppf(1-alpha2, df=f)
    Sw = S*x2w/f
    r1 = r[1]
    # return l, Sl, Sr, Sw, r1
    dataDict = {"t": l.tolist(), "amp": Sl.tolist(), "fre": r1.tolist(), "Sr": Sr.tolist(), "Sw": Sw.tolist()}
    return dataDict